Large-Scale Monte Carlo Study of a Realistic Lattice Model for Gai-^Mn^As. 
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The properties of Mn-doped GaAs are studied at several doping levels and hole compensations, 
using a real-space Hamiltonian on an fee lattice that reproduces the valence bands of undoped GaAs. 
Large-scale Monte Carlo (MC) simulations on a Cray XT3 supercomputer, using up to a thousand 
nodes, were needed to make this effort possible. Our analysis considers both the spin-orbit interac- 
tion and the random distribution of the Mn ions. The hopping amplitudes are functions of the GaAs 
Luttinger parameters. At the coupling J~1.2 eV deduced from photoemission experiments, the MC 
Curie temperature and the shape of the magnetization curves are in agreement with experimental 
results for annealed samples. Although there are sizable differences with mean-field predictions, the 
system is found to be closer to a hole-fluid regime than to localized carriers. 
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PACS numbers: 75.50.Pp, 71.10.Fd, 72.25.Dc 

Introduction. The study of diluted magnetic semi- 
conductors (DMS) is a rapidly expanding area of 
researcbi^i 3 -^ triggered by the discovery of high Curie 
temperatures (Tc) in some of these materials, and by 
their potential role in spintronics devices^ Experimen- 
tally, much progress in the study of DMS was recently 
achieved after annealing techniques reduced the rates 
of compensation, allowing for higher Tc's and for the 
intrinsic properties of DMS materials to be carefully 
investigated^ Theoretical studies have focused on two 
idealized scenarios: (1) Impurity-Band (IB) models; 8 
and (2) Hole-Fluid (HF) models (based on mean-field 
approximations),^ which produce dramatic differences 
in most physical quantities. 

These rather different views (1) and (2) of the same 
problem show that the use of unbiased computational 
techniques, such as lattice MC simulations, is crucial for 
progress in the modeling of DMS materials since it is 
unlikely that the coupling constants and densities locate 
the DMS systems exactly in one extreme or the other. 
Flexible methods that can interpolate between both lim- 
its within the same formalism are needed. However, 
while considerable progress has been achieved in numer- 
ical studies of simplified models^, the complexity of the 
real problem (involving several bands on an fee lattice) 
has prevented its detailed analysis until now. 

Here, we report a comprehensive numerical Monte 
Carlo study of a realistic lattice model for Mn-doped 
GaAs, including spin-orbit coupling, as well as the effects 
of random Mn doping. This large-scale computational ef- 
fort was possible by using the Cray XT3 supercomputer 
operated by the National Center for Computational Sci- 
ences. Our simulations made use of up to 1,000 XT3 
nodes. Parallelization was used in different ways: (i) 
to study different regions of parameter space (densities, 
couplings), and (ii) to average over different configura- 
tions of Mn locations. In all cases, the use of hundreds 
of processors in a single parallel run poses several techni- 



cal challenges that are best handled by supercomputers 
with low latency and scalability, rather than by conven- 
tional clusters of PC's. In fact, this study would have 
taken several years without access to a supercomputer 
with thousands of processors, such as the Cray XT3. 

The Model. We have constructed a real-space fee- 
lattice Hamiltonian whose kinetic-energy term maps into 
the Luttinger-Kohn model^ when fc-space Fourier trans- 
formed and at k — > 0. As a consequence, the hopping am- 
plitudes are functions of (tabulated) Luttinger parame- 
ters, and thus they are precisely known^ To incorporate 
the spin-orbit (SO) interaction, we work in the \j, rrij) ba- 
sis, where j can be 3/2 or 1/2 (since we consider the p 
orbitals, 1=1, relevant at the V point of the GaAs va- 
lence band). Consequently, there are 6 possible values 
for rrij , indicating that this is a fully 6-orbital approach, 
arising from the 3 original p orbitals and the 2 hole spin 
projections. The Hamiltonian is formally given by 
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where a, b take the values =, § (for j=3/2), or | (for 
j=l/2), and a and a' can be 1 or — 1. The Hund term 
describes the interaction between the hole spins si (ex- 
pressed in the \ j, rrij) basis^ 3 .) and the spin of the localized 
Mn ion Si. The latter is considered classical (|Si|=l)ji£ 
since it is large S = 5/2*& /j, + v are the 12 vectors indi- 
cating the 12 nearest-neighbor (NN) sites of each ion lo- 
cated at site i, while I are random sites in the fee lattice. 
Aso=0.341 eV is the spin-orbit interaction strength 17 
in GaAs. The hopping parameters, tjL a , b , are complex 
numbers, whose real and imaginary parts are functions 
of the Luttinger parameters^ 
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FIG. 1: (color online) (a) Curie temperature vs. J, for x=8.5% 
and p~0.75. The MC results are indicated by circles, while 
the continuous line is the MF prediction.—^ Inset: MC re- 
sults for larger values of J to observe the crossover toward 
a localized picture. Vertical lines indicate the experimen- 
tally acceptable range of J. (b) MC calculated Tc vs. p, 
at x—S.5%. The blue dots are experimental results,— and 
the solid line the MF prediction. 



Equation (1) can be studied with standard MC tech- 
niques for systems involving fermions and classical spins. 
These conceptually simple methods were already de- 
scribed in much detail in previous DMS investigations 10 
and also in manganese oxide studies, thus details will not 
be repeated here. Numerically, we analyze clusters that 
contain N — N x N y N z unit cells (JVj is the number of cells 
along the spatial direction i) of side ao (ao=5.64Ajil is 
the GaAs cubic lattice parameter) . Since in an fee lattice 
there are 4 ions associated to each cell, the total number 
of Ga sites is Nc a =AN . Since there are 6 single fermionic 
states per site, the diagonalization of a 6Nc a x 6Nc a ma- 
trix is needed at each step of the MC simulation, which 
demands considerable computational resources for large 
enough clusters. The diagonalization can be performed 
exactly for values of Nc a up to 500. We show below 
that lattices with iVj = 4 (256 sites) are large enough to 
study Mn dopings x and compensations p in the range of 
interest, with sufficient precision for our purposes. Nom- 
inally there should be one hole per Mn ion, but p can be 
smaller than 1 due to hole trapping defects, thus p and 
x are considered independent in this study. 

Results. The highest Tc experimentally observed 
in bulk Gai-^Mn^As is ~ 150 K, at x=8.5% and 
j?~0. 7 19 i 20 . The system is metallic, and the magnetiza- 
tion vs. temperature displays mean- field behavior^ In 
Fig. la, we present the (MC calculated) Tc as a function 
of the coupling J, at x—8.5% and pw0.75. The results 
shown were obtained on lattices containing 256 (Ga,Mn) 
ions, and using an average over at least 5 different disor- 
der configurations (only small differences were observed 
among the Mn configurations). Results for lattices with 



up to 500 sites have also been calculated for some param- 
eters (see below). At the realistic J =1.2 eV—, Fig. la 
shows that the critical temperature is Tc=155 ± 20^. 
Since J is not accurately known, this excellent agreement 
with experiment's, could be partially fortuitous, but at 
least the results indicate that a reasonable quantitative 
estimation of the real Tc can be made via MC simulations 
of lattice models. The solid line in the figure corresponds 
to the MF results The quantitative MF-MC agree- 
ment at small J provides a strong test of the reliability 
of the present MC approach. At J =1.2 eV, the MF T c 
is ~300 K, showing that at these couplings and densities 
appreciable differences between MC and MF exist: the 
fluctuations considered in the MC approach cannot be 
neglected. The inset of Fig. 1 demonstrates that even- 
tually for very large values of J the MF approximation 
breaks down, as expected. The MC simulations show 
that Tc reaches a maximum for J w 12eV, of the order 
of the carriers bandwidth, and then it decreases due to 
the tendency of holes toward strong localization. This 
"up and down" behavior can only be obtained with lat- 
tice MC simulations valid at arbitrary values of JJ£ 

At J~1.2eV the system is closer to a hole-fluid than a 
localized regime as suggested by the magnetization vs. T 
curve, displayed in Fig. 2a. This curve has Curie- Weiss 
shape in qualitative agreement with both experimen- 
tal results' 9 - and previous MF calculations! 9 ' 18 ^ 3 This 
qualitatively correct shape of the magnetization curve 
was not obtained in previous lattice MC simulations. 10 
Size effects are mild as it can be seen in Fig. 3a 
where data for magnetization vs. T are presented for 
x = 8.5%, p « 0.75, and J = 1.2eV in lattices with 
(N x ,N y ,N z ) = (4,3,3), (4,4,4), (5,4,4), and (6,4,4), 
i.e., with iV Ga =144, 256, 320, and 384. Results for 
N=500 were obtained for lower doping (see Fig. 2b). 
Considering together the results for the different size clus- 
ters the estimated Tc «155±15Kis still in agreement 
with experiments (and also with the 256 sites results). 
Regarding the Curie- Weiss (CW) shape of the magneti- 
zation curve, we have phcnomcnologically observed that 
the finite spin-orbit coupling plays a crucial role in this 
respect. In Fig. 3b we show the magnetization vs. T for 
x = 8.5%, p = 0.75, and J = 1.2eV for A so = OMeV 
(squares) and Aso = (circles): only the nonzero SO 
coupling produces CW behavior. In addition, we have 
noticed that the CW shape is also missing with the 
4 orbital (with j = 3 / 2 ) model that results in the limit 
Aso ~~ ¥ oo. 14 This indicates the important role that a 
realistic representation of the valence band plays in prop- 
erly describing the thermodynamic observables. 

The charge distribution in the cluster provides inter- 
esting information. In the HF scenario, the charge is 
assumed to be uniformly distributed while in the IB pic- 
ture the charge is strongly localized. Fig. 3c indicates 
that in the realistic regime with J = 1.2 eV, x = 8.5%, 
and p w 0.75 the charge is fairly uniformly distributed. 
The slightly darker points correspond to the sites where 
the Mn are located. They have charge of the order of 
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FIG. 2: (color online) (a) Magnetization M vs. T, for £=8.5% 
and several p's (indicated), using a 256 sites lattice (open 
symbols). Averages over 5 Mn-disorder configurations are 
shown, (b) Same as (a), but for x=3%. Close circles are 
results for a 500 sites la ttice. The magnetization is mea- 
sured as M — VM ■ M, with M the vectorial magnetiza- 
tion. As a consequence, for fully disordered spins, A4 is still 



at large temperatures (M(T — > oo) = l/y/xNoa) unrelated 
to ferromagnetism. Thus, we plotted M = (M — M(T — » 
oo))/(l — M(T — » oo)), i.e. the background was substracted. 



20% above the MF value defined as n.MF=n/Nc a (with 
n the number of holes) . As shown in Fig. 3d, charge lo- 
calization occurs when J is increased to large values such 
as 12 eV. The dark circles at the Mn sites have charge 
intensities about 20 times the MF value, with very little 
charge found away from the impurities. 

The lack of charge localization effects in the ordered 
state is concomitant with the absence of a notorious im- 
purity band in the density-of-states (DOS) (Fig. 4a). In- 
creasing J, an IB regime is eventually observed, given 
confidence that the study is truly unbiased. This occurs 
for J « 4eV and beyond, with the IB becoming totally 
detached from the valence band at J « 16eV. Figures 3c, 
d and 4a show that the degree of spatial hole localization 
is correlated with the development of the IB. In addition, 
when the holes are localized, the M vs. T curves present 
substantial deviations from MF behavior (not shown), 
with a different concavity as that of Fig. 2aJ£ 

It is well known that the carrier density in DMS is 
strongly dependent on sample preparation. Due to de- 
fects, p is much smaller than 1 in most samples. In 
Fig. lb, the MC calculated T c vs p for £=8.5% and 
J=1.2ey is shown. Tq increases with p, and it reaches a 
maximum at p~l as in previous theoretical^^ and ex- 
perimental results^ ' 19 ! 25 ! 27 also shown in the figure. The 
agreement between MC results and experiments is once 
again quite satisfactory. The figure suggests that if p=l 
were reachable experimentally, then Tq could be as high 
as ^200 K. We also observed qualitative changes in the 



FIG. 3: (color online) (a) M (defined as in Fig. 2) vs T for 
different lattice sizes for x = 8.5%, p w 0.75, and J = 1.2eV; 
(b) M vs. T for the same parameters as in (a) on a 256 sites 
lattice with (without) spin-orbit interaction indicated by the 
squares (circles); (c) Charge density normalized to the MF 
value (see text), for x = 8.5%, pw 0.75, T=10K, on a 256- 
sites cluster for J—1.2eV. The color intensity is proportional 
to the charge density (see scale), (d) Same as (c), but for 
J=12eV. 
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FIG. 4: (color online) (a) Density-of-states, for £=8.5%, 
p«0.75, and several J's. The dashed vertical lines indicate 
the position of the chemical potentials, (b) Same as (a) but 
for 1 = 3%. 



magnetization curve varying p (Fig. 2a): as p is reduced, 
the magnetization changes from a Brillouin form to an 
approximately linear shape with T . In fact, the observed 
p dependence is once again similar as found experimen- 
tally, with p modified by annealing^ the Mn disorder 
plays a more dominant role when the number of holes is 
reduced. In spite of the deviations from MF behavior at 
small p, we could not observe a clear IB as p was reduced, 
for the coupling J used in our analysis. 
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Consider now the low Mn-doping regime. The latest 
experimental results suggest that (Mn,Ga)As should still 
be metallic at x=3%. 25 A metal-insulator transition is 
expected at x < but MC studies for such very 

small dopings need much larger clusters than currently 
possible. Our MC simulations indicate that at low Mn- 
doping, the dependence of Tq with J is similar as in the 
higher doping case. 

No clearly formed IB is observed in the x=3% DOS 
displayed in Fig. 4b: the IB's are formed for similar val- 
ues of J with increasing J, at both x=3% and 8.5%, 
as deduced from an analysis of fermionic eigenvalues for 
spin-ordered configurations! 14 i 29 

Summary. Here, a MC study of an fee lattice model for 
DMS compounds, including the realistic valence bands of 
GaAs, the spin-orbit interaction, and the random distri- 
bution of Mn dopants has been presented. The use of the 
Cray XT3 at ORNL made this effort possible. The results 
show magnetizations and Ic's in reasonable agreement 
with experiments. The simulations show that the carri- 



ers tend to spread over the entire lattice, and they reside 
in the valence band at realistic couplings, qualitatively in 
agreement with MF— £ and first principles 5 calculations 
in the same parameter regime, as well as with experimen- 
tal data on annealed samples^^ However, an IB band 
populated by a fraction 1 — p of trapped holes that do not 
participate in the transport properties is not ruled out by 
our results. The MC method described here opens a new 
semi-quantitative window for theoretical research on the 
properties of DMS materials. 

Acknowledgments. We acknowledge discussions with 
T. Dietl, A. MacDonald, J. Sinova, T. Schulthess, F. 
Popescu and J. Moreno. Y.Y., A.M., and E.D. are sup- 
ported by NSF under grants DMR-0443144 and DMR- 
0454504. G.A. is sponsored by the Division of Materi- 
als Sciences and Engineering, Office of Basic Energy Sci- 
ences, U.S. Department of Energy, under Engineering, 
Office of Basic Energy Sciences, DOE, under contract 
DE-AC05-00OR22725 with ORNL, managed and oper- 
ated by UT-Battelle, LLC. 



1 H. Ohno, Science 281, 951 (1998). 

2 T. Dietl, et al, Science, 287, 1019 (2000). 

3 T. Jungwirth et al, Rev. Mod. Phys. 78, 809 (2006). 

4 A. MacDonald et al, Nature Materials 4, 195 (2005). 

5 T. Schulthess et al, Nature Materials 4, 838 (2005). 

6 I. Zutic et al, Rev. Mod. Phys. 76, 323 (2004). 

7 S. J. Potashnik et al, Appl. Phys. Lett. 79, 1495 (2001). 

8 M. Berciu and R. Bhatt, Phys. Rev. Lett. 87, 107203 
(2001). 

9 M. Abolfath, T. Jungwirth, J. Brum and A. MacDonald, 
Phys. Rev. B63 054418 (2001). 

10 One-band models were studied in G. Alvarez et al, Phys. 
Rev. Lett. 89, 277202 (2002) and M. Mayr et al, Phys. Rev. 
B65, 241202(R) (2002); two-band models in F. Popescu 
et al, Phys.Rev.B73, 075206 (2006); MC results for the 
kinetic-exchange model are in J. Schliemann et al, Phys. 
Rev. B64, 165201 (2001); and LSDA-based estimations of 
T c are in J. Xu et al, Phys. Rev. Lett. 84, 097201 (2005). 

11 J.M. Luttinger and W. Kohn, Phys. Rev. 97, 869 (1955). 

12 J.M. Luttinger, Phys. Rev. 102, 1030 (1956). The hoppings 
used in our study are all of order 1 eV. 

13 These are 6 x 6 matrices; see Eq.(A.lO) (in Ref@) andQJ. 

14 Details will be published elsewhere (in preparation). 
The classical limit corresponds to limh_o,s , ^oo frS = fio = 
6.58 x 10 -16 eV s. Since the parameter J, experimentally 
measured, is proportional to h via the Bohr magneton, 
it results that J = Kh where K is a constant. Thus, 
linifc^o.s^oo JS = limft^o.s^oo KhS = Kh = J. See G. 
Q. Pellegrino, K. Furuya, and M.C. Nemes, Chaos 5 463 
(1995); G.Q. Pellegrino, K. Furuya, and M.C. Nemes, Re- 
vista Brasileira de Ensino de Fisica 20, 321 (1998); L.G. 
Yaffe, Rev. Mod. Phys. 54, 407 (1982). This approach is 
also supported by comparisons of our numerical data with 
MF calculations for the quantum case. 

16 J. Szczytko et al, Phys. Rev.B60, 8304 (1999); M. Lin- 
narsson et al, Phys. Rev.B55, 6938 (1997); O.M. Fedorych 



et al, Phys. Rev.B66, 045201 (2002). 

17 Peter Yu and Manuel Cardona, Fundamentals of Semicon- 
ductors, Fhird Edition, Springer- Verlag. 

18 T. Dietl, private communication. 

19 K.C. Ku et al, Appl. Phys. Lett. 82, 2302 (2003). 

20 T. Jungwirth et al, Phys.Rev.B72, 165204 (2005); KY. 
Wang et al, AIP Conf. Proc.772, 333 (2005). 

21 N=500 results were obtained for lower doping (Fig. 2b). 

22 J. Okabayashi et al, Phys. Rev. B58, R4211 (1998). 

23 T. Dietl et al, Phys.Rev.B63, 195205 (2001). 

24 John Schliemann et al, Appl. Phys. Lett. 78, 1550 (2001). 

25 KW. Edmonds et al, Appl. Phys. Lett. 81, 3010 (2002). 

26 S.-R. Yang and A.H. MacDonald, Phys. Rev. B67, 155202 
(2003). 

27 S.J. Potashnik et al, Phys. Rev. B66, 012408 (2002). 

28 J. Okabayashi et al, Phys. Rev. B64, 125304 (2001); J. 
Okabayashi et al, Physica E10, 192 (2001). 

29 Coulomb interaction. While the J term in Eq.(l) is crucial 
for ferromagnetism, it is expected that most of the carrier 
binding energy, at least at very low doping, originates from 
the attraction between the Mn ions and the hole carriers. 
To mimic this attractive Coulomb interaction, we added to 
Eq.(l) a term of the form Hcoui = F^ t ni, using values 
of V up to 2 eV, which is needed to generate the ~0.1 eV 
binding energy of (Mn,Ga)As. At small J, V tends to ef- 
fectively increase J, since it increases the chances for the 
carriers to interact locally with the Mn spins. As a conse- 
quence, adding a Coulomb attraction as described above, 
does not generate a case at £=8.5% or 3%, presenting both 
an IB and a Tc close to the experimental range. It is pos- 
sible though that the Coulomb attraction develops an IB 
for smaller values of x (in the insulating regime^) that 
cannot be studied with our finite clusters. Details will be 
presented elsewhere.— 

30 F. Popescu et al, in preparation. 



